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Abstract 

The standard approach for integrating the multidimensional Vlasov equation using 
grid based, conservative schemes is based on a time splitting approach. Here, we 
show that although the truncation error is of second order, time splitting can intro- 
duce systematic heating of the plasma. We introduce a backsubstitution method, 
which not only avoids this deficiency but also is computationally less expensive. The 
general approach is demonstrated in conjunction with Boris' scheme for evaluating 
the characteristics. 
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1 Introduction 

Vlasov's equation is fundamental for numerous problems in plasma theory. 
This kinetic equation describes the behaviour of the single particle distribu- 
tion functions of a collisionless plasma under the influence of electric and mag- 
netic fields. Coupled with the equations for the electromagnetic fields and the 
evaluation of the moments of the distribution functions one obtains a highly 
nonlinear system of differential and integral equations. Only a few very sim- 
ple problems can be solved analytically. For this reason numerical simulations 
of Vlasov's equation have become an important tool for theoretical plasma 
physics. 

One type of computer simulation approach integrates the distribution func- 
tion directly on a high-dimensional numerical grid in phase space. Here one 
dimension is needed for every space component and for every velocity compo- 
nent. Following the original work by Cheng and Knorr [1] much progress has 
been made on improving the accuracy of the advection schemes. 
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The integration of the distribution function can be carried out in a number 
of different ways: The simplest schemes are finite difference schemes. They 
are relatively easy to implement but suffer from numerical instabilities and 
nonpositivity of the distribution function. Conservation laws such as the con- 
servation of particle number can be implemented but complicate the scheme 
greatly [2]. Conservative methods, on the other hand, discretise the distribu- 
tion function by integrating over the numerical grid cells [3]. The advantage of 
these methods lies in the fact that the particle number is naturally conserved 
and no artificial sources or sinks of particles are introduced. Semi-Lagrangian 
methods (e.g. [4]) follow the characteristics backwards and interpolate the dis- 
tribution function at the origin of the characteristic. The interpolated value 
is then transported forward to the grid points. Semi-Lagrangian methods do 
not naturally conserve the particle number but can easily be made to preserve 
positivity. 

Most of the above methods are, however, developed for a one-dimensional 
advection problem. When used for a one-dimensional electrostatic system in 
which the physical phase space is two-dimensional, a time splitting method 
is employed which was already proposed in [1]. Although Semi-Lagrangian 
methods in principle allow to integrate the distribution function directly on 
the high-dimensional grid, the time splitting technique is also used to simplify 
the computation [5]. The general idea is that in higher dimensions this time 
splitting can be generalised in a straightforward way [6] . 

We will show in this paper that, when including a magnetic field, this simple 
time splitting — although second order — can cause dissipation due to errors 
which are always in the same direction. This implies that the temperature 
of the system will increase systematically. We will also present an alterna- 
tive method which we named backsubstitution method. The backsubstitution 
method not only eliminates the problems of the time splitting method but is 
also computationally less expensive. 

In section 2 we will present the basic underlying equations. Section 3 will 
describe the time splitting method and show how a systematic error develops. 
In section 4 we will present the backsubstitution method which we will apply 
to Boris' scheme in section 5. Section 6 discusses simulations of Bernstein 
waves using the different schemes to provide a comparison. Section 7 gives 
some concluding remarks. 



2 General Problem 

The basis of the kinetic plasma description is the distribution func- 
tion /(x, v,i) which expresses the particle density in phase space. Here 
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/(x, v.t) d^x d^v is the number of particles in a phase space volume d^x d^v 
located at (x, v) at time t. In a collisionless plasma the evolution of the dis- 
tribution function is given by Vlasov's equation 



^ + V • V./ + ^ (E + V X B) • Vv/ = 0, (1) 
at m 

where E and B are the electric and magnetic fields which have to be deter- 
mined self-consistently. Vlasov's equation describes the advection of values of 
the distribution function along particle characteristics given by Newton's law 
of motion. 

One central property of Vlasov's equation is the conservation of the phase 
space density, which directly translates into a conservation of mass and charge 
in a closed system. For this reason it is natural to use a conservative scheme for 
simulating Vlasov's equation (for 1-dimensional schemes, see e.g. [7]). Today, a 
diversity of Eulerian schemes, all with high accuracy and different advantages 
and disadvantages, are available (see e.g. [5,8] and references therein). These 
schemes normally solve the one-dimensional advection problem, 

dtf{x,t)+udj{x,t) = 0. (2) 

By integrating over a finite time step one obtains 

I f{xX^')dx^ I f{xX)dx. (3) 

Here X{s,t,^) denotes the characteristic with parameter s that satisfies 

For the one dimensional, electrostatic Vlasov-Problem 

% + vdJ + ^EdJ = 0, (4) 
ot m 

a splitting technique is then usually employed. Here one integrates the advec- 
tion in the x-direction by A.t/2, then in i>-direction by At and then again 
in x-direction by At/2. This produces a second order scheme which can be 
written as 

T^{At/2)T^{At)T^{At/2). (5) 
Here denotes the numeric advection operator in the /c-dimension. 
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3 Time splitting 



The success of the time sphtting for the one dimensional electrostatic prob- 
lem motivates a common suggestion to extend the splitting technique to treat 
higher-dimensional systems. Since the spatial dimensions are completely in- 
dependent of each other this results in the following second order scheme for 
the full three dimensional system 



T,(At/2)T,(At/2)r,(At/2) 

T,, ( At/4)r,^ ( At/2)T,, ( At/4)T,^ ( At)T„, ( At/4)T,^ {M/2)T,^ (At/4) 



In each of these sub-steps a one dimensional transport equation of type (2) 
is solved. For each of these equations the characteristics are calculated and 
then projected onto the corresponding direction. This implies that even for 
a hypothetical exact one-dimensional integration scheme the characteristics 
are still only approximated by a second order time splitting scheme. For the 
following discussion, we will consider only the velocity part of the integration 
scheme since this determines how well the particle temperatures are described. 

For a purely electrostatic system the above integration behaves well and errors 
only occur due to repeated application of the advection scheme. The reason for 
this is the independence of the change of the velocity component Af ^ on the 
velocity Vk- With a magnetic field, however, the change of velocity Af^ over a 
finite time step does depend on the velocity Vk- Here we assume that the inte- 
gration scheme for the characteristics is at least second order. To investigate 
the error caused by this method, we take the exact characteristic in v-space 
and approximate it using the time splitting scheme. During integration of the 
characteristic, the electromagnetic fields are assumed to be constant. Without 
loss of generality we let B = Bt. and move the origin in velocity space to 
Vo = E X B. For simplicity we assume — 0. Ez would only add a 
constant acceleration in the t'^-direction, and leads to the same result. 

In this setup the characteristics in velocity space are simple concentric circles 
around the origin and we can neglect the Vz coordinate completely. During 
the time interval At the whole Vx,Vy-plane rotates by an angle 4> = AtqB/m. 
Taking a velocity 



we split the rotation into three steps according to the time-splitting scheme 



T,{At/2)Ty{At/2)Tz{At/2). (6) 



V = {V:„Vy) 



(7) 



T,,(At/2)T,^(At)T„,(At/2). 
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This results in the following 



v" = {v^ cos(0/2) - Vy sin(0/2), Vy) , (9) 
v''=«,<cos(0)+<sin(0)), (10) 
= (r;^os(0/2) - t;^ sin(0/2), t;^) . (11) 

Inserting v" into and then into v'^^^™ results in a lengthy expression for 
^iicw Xaking the norm of v'^'^™ and expanding this expression for small angles 
(f), i.e. small time steps At, gives 

{v-f-vl + vl-{vl + iy' + 0{<l>% (12) 

By construction this is, of course, second order in 0. However, one can sec that 
the second order error is always negative and thus introduces a systematic 
error. 

The time splitting method (8) can also be interpreted as performing the indi- 
vidual steps (9)-(ll) in first order in 0. This corresponds to 

y''^{V,-Vy(l>/2,Vy) , (13) 

v^=«,< + <0), (14) 
v-- = (^^-^j0/2,^J). (15) 

With respect to equation (2) this scheme is obtained by holding u constant 
for each step. Taking the square of this v"^^^ results in 

v-^-'^vl + vl + ^v^Vy + Oi<f>% (16) 

One can observe that in this case the second order disappears, and the third 
order error is not systematic, but depends on the signs of and Vy. In the 
following we will refer to eqs (9)-(ll) as scheme A and (13)-(15) as scheme 
B. 

In Fig. 1 the error of the magnitude | v| of the velocity after a quarter gyration 
is plotted against the rotation angle Aip of the individual step. The error is 
normalised to the initial velocity. The solid hne shows the result of scheme A 
while the dashed line shows the result of scheme B. The dotted line represents 
the result of scheme B with alternating order of the v^-Vy integration. For 
scheme A, a total error of about 2.5% is accumulated after a quarter gyration 
when Aip 0.045. After a full gyration the error sums up to 10% (not shown). 
This value of Acp corresponds to roughly 140 integration steps for the full 
circle. Using less steps, i.e. larger A(f results in even larger errors. 

To understand the direction of the error we note, that in eq. (3) the charac- 
teristics are integrated backwards from a grid point to obtain the source of 
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Fig. 1. Error Verr of the magnitude of the velocity after a quarter gyration depending 
on the angle of an individual step. Curves are plotted for the schemes A, B and 
scheme B with alternating order of the Vx-^y integration. 

the distribution function for that grid point. The distribution function is then 
transported from that source to the grid in some manner that depends on the 
numerical scheme. The negative sign in the second order of eq. (12) imphes 
that the source is always located closer to the rotation centre than the grid 
point. Thus the values of the distribution function arc transported outwards 
from the rotation centre. This results in an effective heating of the distribution 
function. 

Using scheme B the errors are smaller but not zero. Here an error of 2% 
is observed when Aip k, 0.4. This is equivalent to roughly 16 steps for a 
full gyration. Because the direction of the error in scheme B depends on the 
values of Vx and Vy, one can further increase the accuracy by alternating the 
order of the splitting. In the two dimensional case considered here this simply 
implies alternating the roles of Vx and Vy. Using the alternating scheme the 
overall error in the velocity magnitude is reduced to almost zero. However, 
when looking at the relative phase error after a quarter gyration (Fig. 2) 
no significant improvement can be observed. While scheme A still shows the 
largest error, the errors for scheme B with and without alternating oder of 
integration are roughly comparable up to a A</7 of 0.5. For this value of Ac/? 
the phase error is approximately 1%. 

When the Vlasov equation is solved on a discretised grid errors are worse 
but the main sources of these errors are highlighted by the above analytical 
argument. 
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Fig. 2. Error (perr of the phase of the velocity after a quarter gyration depending on 
the angle of an individual step. Curves are plotted for the schemes A, B and scheme 
B with alternating order of the v^-Vy integration. 

4 Backsubstitution 

In this section we want to present an alternative method for integrating 
Vlasov's equation that does not suffer from the above drawbacks. Here we 
will present first the general idea of this backsubstitution method and then 
write down the equations for the general system described above. 

Suppose we are given a one dimensional integration scheme for the transport 
eq. (2). To create a scheme for the integration of the three-dimensional velocity 
space there is no other choice but to split the full three-dimensional problem 
into a number of one-dimensional substeps. For each of these substeps the 
characteristics will be calculated and then projected onto the direction of the 
advection step. We still have the freedom, which characteristics to integrate 
and in which order to integrate them. 

To start with, let us again consider the standard case described in the last 
section. Our aim is to formulate a splitting scheme in which the characteristics 
are integrated exactly, and which uses the minimum number of integration 
steps. Since we can ignore the t'^-direction, this means we want only two 
integration steps, one for v^, and one for Vy. 

The distribution function is first shifted in the v.^, and then in the Vy direction. 
Figure 3 illustrates the first step while Figure 4 illustrates the second step. 
Both shifts together should transport the value of the distribution function 
from a source point S — {S^, Sy) of a characteristic to its destination point 
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Fig. 3. Schematic diagram of the first step of integration of the characteristics 
using the backsubstitution algorithm. The distribution function is shifted in the 

i;2:-direction from the source point 5'-^) of the characteristic to the grid point G. 
Gray hues show corresponding characteristics in the other cehs. The dashed hne is 
the characteristic ending on G which is important in the second step. 

D = {Dxi Dy) with S = X(t — At, t, D). Here the indices and z are used 
to denote the velocity components v^-, Vy and Vz- This means that we aim to 
find a scheme such that 

r^{D,,Dy)^r'\S,,Sy). (17) 

In the first step the shift in has to transport / from S to an intermediate 
point {Dx^ Sy). In the semi-Lagrangian schemes which we are considering here, 
the characteristics are integrated backward from the grid points. This implies 
that in the first step (1) the grid point G has to coincide with the intermediate 
point G = {D'i\S^^^), or equivalently S^^^ = {S^^\ Gy) and L>(i) = (G^, 
We have displayed these characteristics in Figure 3. In this way the distribution 
function has been shifted along Vx according to 

r'^\Gx.Gy)^r'\S'^^\Gy). (18) 

Given a sufficiently smooth behaviour of the characteristics we can assume that 
the interpolation scheme causes all other points of the distribution function 
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Fig. 4. Schematic diagram of the second step of integration of the characteristics 
using the backsubstitution algorithm. The distribution function is shifted in the 

direction from the intermediate point to the destination point G of the charac- 
teristic. Gray hnes show corresponding characteristics in the other cells. 

to be shifted accordingly. This is particularly true for the characteristic that 
ends in the grid point G (dashed line in Figure 3). This characteristic will be 
important in the following step. 

In the second step (Figure 4) we, therefore, need to choose the characteristic 
that ends in G. Then the source point S**^^-* is given by S^"^^ = X{t — At, t, G). 
The shift is performed in the Vy direction so that 

r^{Gx,Gy)^r'^^{Gx,Sf^). (19) 

Since in the first step we had (assuming again correct interpolation) 

r*-(G.,5f) = r<i(5(^),5f), (20) 

we finally have 

r^{Gx,Gy)^r'\sf\s^^'^). (21) 

We now use this motivation to write down a general scheme for three- 
dimensional velocity space. For every grid point G we perform the integra- 
tion in three one-dimensional substeps, one for each component Vx,Vy,Vz- For 
each integration a source coordinate S^^\S^^^ and S^^^ is calculated from a 
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characteristic which docs not necessarily pass through G. To find S^"^ for the 
T;j;-integration we demand 

Gx = D^}) , (22) 
(x(t-At,t,D«))^ , (23) 

G,= {x{t-At,t,D^'y))^ . (24) 

In general this is a nonhnear system of equations for the components D^^^ 
and D^^^ . The details of this system depend on the way the characteristics are 
calculated. Given D^^^ one then has 

5« = (x(t-At,t,D«))^ , (25) 

and the integration can be performed from to G^ in the i'j;-direction. 
Similarly we demand for the v^^-integration 

G. = , (26) 
Gy = , (27) 
G,= {x{t-At,t,D('^))^ . (28) 

D^"^^ differs from G only in the t'^-component. Once D^^ is found we have 

5(2)= (^X(t-At,t,D(2)))^. (29) 

Again the integration is now performed from S^^^ to Gy in the v^-direction. 

The T;2-integration finally is straightforward. Since D^^^ — G we have 

Si'^ = {X{t-At,t,G))^ (30) 

and the integration is performed from S^^^ to G^ in the T;2-direction. 

We want to emphasise that this scheme integrates the characteristics exactly, 
which means that in terms of Figs. 1 and 2 the backsubstitution scheme has 
errors which are exactly zero. 

5 Application to Boris scheme 



The main task now is to calculate the characteristics or their approximations 
in the presence of a magnetic field. A commonly used approach is the Boris 
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scheme [9]. Here the integration step is formulated as an implicit finite differ- 
ence scheme 



V 



n+l 



At 



n+l 



X B 



m 



The electric and magnetic forces are separated, 

2 m 



V ' = V 



n+l 



(31) 

(32) 
(33) 



leading to 

v+ — v~ q 
At ~ 2m 

The transformation from v~ to v+ is a pure rotation with an angle 9 where 



(v+ + v") X B 



(34) 



tan ■ 
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At qB 



(35) 



For small angles 9 this is close to the exact angle ^exact = AtqB/m. 



In contrast to the original scheme of Boris, we aim to trace the characteristics 
backward in time. This means we want to find v~ in terms of v"^. We thus 
reverse the original scheme and rotate v"*" by To implement this rotation 
the vectors t and s are defined 



t = 



-At 



2t 



s = 



(36) 



2m' 1 + ^2- 

Then the rotation is performed in two steps 

v' = v+ + v+ X t , (37) 

and 

V" = v+ + v' X s . (38) 

This scheme now supphes -S" = v" = i'^x^'^y^'^z) terms of D = v"+^ = 
(f""*"^, fy"*"^, i'""'"^). To facihtate the further calculations we insert (37) into 
(38) and separate v~ into it's components 

V~ ^ {1- Syty - S:,tz)v^ + {Sytx + S z)Vy + {Szt^ 



{Sxty - S^)V^ + (1 - S^tx - S^tz)v^ + {S^ty + S^)v^ 
{Sxtz + Sy)v^ + {Syt^ - Sx)v:^ + (1 - Sxtx - Syty)vi 



(39) 
(40) 
(41) 



We now need to solve the systems of equations (22 - 24) and (26 - 28) for the 
first and the second backsubstitution step. As stated before, the third step 
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is straightforward since D^^^ is already known. The complete problem can be 
written in the form 



nf n+l 



X \ ^X 1 



<, , (42) 

< = <K+\<+\^':), (43) 

< = <«+\<^\<+')- (44) 

Since the bijections between and f ~, on one hand, and and f +, on the 

other hand, are trivial (see eqs (32) and (33)) it is sufficient to formulate the 
three steps 

% ^'^xi'^t^v-^^7) ^ (45) 

Vy ^v-{vt,v^,v-) , (46) 

^'^ziyt^v^.vt) . (47) 



To find (45) we take eqs (40) and (41) and solve for and v'^ giving 

(1 - Sxtx - Syty)V~ - {Sx + Szty)V~ + {sjy{m " l) " S^Hy + S^Sy + S z)V + 

V — 

^ m{sxtx - 1) + 1 + Sx{sx -t^- n^) 

(48) 

+ _ (gg; - Syt^)v- + (1 - - s^t^)v- + {sj^{m - 1) - g^n^ -Sy + SxS^)v+ 

^ m{Sxtx - 1) + 1 + Sx{Sx -U- Tlx) 

(49) 

where m = s ■ t and n = s x t. These can be inserted into (39) which then 
provides the expression (45) for the first step. 

For the second step (46) only eq (41) has to be solved for v+ giving 

(50) 



+ _ -{Sy + SxQv^ + {Sx - Syt,)v+ + VI 



1 ^x^x ^y^y 

With this, vj" can be substituted in eq (40) giving v~ in the form (46). 

By virtue of equation (41) the z-component is already given in the form (47). 
Thus, no further calculation has to be done for the third step. 

Finally we want to discuss the error of Boris' scheme combined with the back- 
substitution method. We again investigate the same problem as formulated in 
section 3 where the velocity vector is rotated around the origin. While Boris' 
scheme introduces a phase error in this rotation, the magnitude of the velocity 
is preserved. Using this combined scheme in a grid based Vlasov solver implies 
that the only diffusion in the system originates from the reconstruction of the 
distribution function. 
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6 Bernstein Waves 



We have applied the schemes described above to the simulation of Bernstein 
waves in a periodic system. These are electrostatic waves propagating at a right 
angle to a given constant magnetic field [10,11]. The ions are treated as a static 
neutralising background, while the electrons oscillate in the electrostatic field. 
We assume that B = B^z and the wavevector k = k^^. Then the dispersion 
relation can be written as 

with 

[n'^ — w^) {[n + 1) —w^] 

Here we used w = u/Q, k = Vthk^/^, = eBz/uif, is the electron cyclotron 
frequency and uip^, is the electron plasma frequency, mg is the electron mass 
and e is the electron charge. F is the gamma function and M is Rummer's 
confluent hypergeometric function. We chose cj^g = Vt^ for all simulations. For 
a given /c^,, the above dispersion relation has an infinite number of solutions for 
uj. We performed the simulations in one space and three velocity dimensions, 
{x^VxiVy^Vz). Although two velocity dimensions would be sufficient for this 
system, we keep the i^^-dimension to make the results transferable to electro- 
magnetic simulations in which the magnetic field is not fixed. The simulation 
box has a length L which was resolved with 64 grid cells. The velocity space 
was sampled with 50 grid cells in each direction in the interval from —-ivth 
to +^Vth ■ The length of the box is chosen to L = 27i/kx so that exactly one 
wavelength of the Bernstein mode fits into the system. In this way the size 
of a grid cell in space is Ax = 27v/64kx. The timestep was chosen such that 
the CFL-condition is satisfied At = Ax/5vth- For the integration of the dis- 
tribution function on the grid we use a fiux conservative and positive scheme 
[7]. 

The first simulation was initialised with the Bernstein wave of the lowest 
frequency mode 1 < uj/Q < 2. Runs were performed for different values of kx 
and using the different integration schemes. The frequency of the wave was 
then determined using a Fourier analysis. The results are shown in Fig. 5 for 
the time-splitting scheme A, time-splitting scheme B, backsubstitution using 
the Boris scheme and backsubstitution using the exact characteristics. We 
can observe that the time-splitting scheme A clearly shows the largest error 
in the dispersion of the waves. The errors of all the other schemes appear 
comparable and are very good for all values oi k > 0.15. The larger errors for 
smaller wavenumbers are due to the choice of the timestep At. Inserting the 
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Fig. 5. Comparison of the dispersion of Bernstein waves between the exact result 
and the Vlasov simulation using different schemes for integrating the characteristics. 
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Fig. 6. Comparison of the dispersion of Bernstein waves between the exact result 
and the Vlasov simulation using different schemes for integrating the characteristics. 

definitions for Ax and k one finds that 



At = 27rO-^ 



320A; 



(53) 



For k = 0.05 this means At = 2Trfl ^/16 or 16 steps for one gyration. 



In another simulation run, the second Bernstein mode 2 < uj/Q < 3 was 
initialised. Fig. 6 shows the result of the time-splitting scheme B and the 
backsubstitution method with exact integration of the characteristics. The 
dispersion relation is shown for values of A; < 0.2, where the errors are largest 
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Scheme 


Comp. time / min 


OJJllLLlIlg OClitJIlitJ I\ 


1 Q9 


Splitting Scheme B 


165 


Backsubs. Boris 


94 


Backsubs. Exact 


137 



Table 1 

Computational time used for a typical run 

due to the choice of the timestep. For this case we observe that the backsub- 
stitution show superior results when compared to the time-splitting method. 

Finally, we want to look at the computational time used by the different 
schemes. Table 1 shows the times used for a typical run. The runs for the dif- 
ferent schemes were carried out with exactly the same conditions on the same 
machine. Here we find a clear advantage of the backsubstitution scheme over 
the sphtting schemes. The backsubstitution scheme with Boris integration of 
the characteristics reduces the computational effort by more than 50% when 
compared to the time splitting scheme A. With time-splitting scheme B this 
improvement is still approximately 43%. The timing for the exact backsub- 
stitution shows less improvement due to the fact that trigonometric functions 
have to be evaluated. The reason for the speed-up is the fact that the backsub- 
stitution method has to integrate the distribution function only once for each 
velocity dimension v^, Vy and Vz- The splitting schemes, on the other hand, 
have to integrate the distribution function 7 times. Although the numerical 
effort of integrating the characteristics in each step is considerably smaller in 
the splitting scheme B, this is only a part of the computational time spent. 
Other parts involve the interpolation of the distribution function and the cal- 
culation of fluxes across the cell boundaries. Considering all the above results, 
the backsubstitution method together with the Boris scheme can be taken as a 
good alternative to the traditional time splitting method if speed is the major 
issue. To obtain the most accurate results, the backsubstitution method to- 
gether with the exact integration of the characteristics is the superior scheme. 
In addition it also is shghtly faster than the time-sphtting scheme. 



7 Conclusions 

We have shown that the time splitting method for integrating Vlasov's equa- 
tion in higher dimensions can introduce systematic errors when used in the 
presence of a magnetic field. These errors originate from the effective split- 
ting of the integration of the characteristics, when a higher order integration 
scheme is used. The errors cause the temperature of the distribution function 
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to increase over time, and thus artificially feed energy into the system. 

The backsubstitution method presented here for the general case of arbitrary 
integration schemes of the characteristics eliminates this problem. Here not 
those characteristics that pass through the grid point are integrated, but those 
characteristics that will give a consistent scheme when executed in sequence 
for the full timestep. This not only provides the best accuracy possible but 
also reduces the number of integration steps. While in three dimensional ve- 
locity space the time-splitting scheme consists of 7 steps, the backsubstitu- 
tion method only uses 3 steps since each component needs to be integrated 
only once. Due to this advantage the backsubstitution method together with 
the Boris scheme typically decreases the computational effort by over 40% as 
compared to a simple time splitting method while the errors remain small. On 
the other hand, highest accuracy can be achieved with the backsubstitution 
method together with the exact integration of the characteristics. 
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